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ABSTRACT 

Although accretion onto supermassive black holes in other galaxies is seen to produce powerful jets 
in X-ray and radio, no convincing detection has ever been made of a kpc-scale jet in the Milky Way. 
The recently discovered pair of 10 kpc tall gamma-ray bubbles in our Galaxy may be signs of earlier 
jet activity from the central black hole. In this paper, we identify a gamma-ray cocoon feature in the 
southern bubble, a jet-like feature along the cocoon's axis of symmetry, and another directly opposite 
the Galactic center in the north. Both the cocoon and jet-like feature have a hard spectrum with 
spectral index ~ —2 from 1 to 100 GeV, with a cocoon total luminosity of (5.5 ± 0.45) x 10 35 and 
luminosity of the jet-like feature of (1.8 ± 0.35) x 10 35 erg/s at 1 — 100 GeV. If confirmed, these jets 
are the first resolved gamma-ray jets ever seen. 

Subject headings: galaxies: active — galaxies: starburst — gamma rays — ISM: jets and outflows 



1. INTRODUCTION 

Supermassive black holes (SMBHs) of 10 6 to 10 10 so- 
lar masses are believed to lie at the center of most 
galaxies and are fed by accretion of ambient gas and 
stars. Accretion-powered jets have been observed at 
various astron omical scales ranging f rom active galactic 
nuclei (AGN) ([Bridle fc Perlev 1984 ). especially blazars 
(BL Lac objects and flat spectrum radio quasars) at 
the bright end, to gamma-ray bursts and Galactic bi- 
naries (stellar mass black holes, neutron stars and cat- 
aclysmic variables). Some of these objects appear to 
have produced jets nearly continuously for at least tens 
of millions years. The mechanism by which jets turn 
on and off is one of the major puzzles in high en- 
ergy ast rophysics, and may be connect ed to star for- 
mation (jAntonuccio-Delogu fe Silkl 120081 ). The relativis- 
ts jets inject significant amounts of energy into the 
medium within which they propagate, creating an ex- 
tended, under-dense and hot cocoon. After decades of 
study, we still lack a complete understanding of the main 
mechanism launching, accelerating, and collimating jets, 
with limited knowledge of the energy content, the com- 
position, and the particle acceleration mechanisms of 
the j ets ( Blandford fc Znajekl 1977 : Blandford fc Pavne 
1982). 



The SMBH at the center of the Milky Way (MW) is 
surrounded by clusters of young star s and giant molec- 
ular clouds ( Morris fc Serabvnl 19961). Although ther e 



are indications of past activity (punvaev et al. 

the SMBH is currently in a quiescent state. Despite 
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the abundant observational evidence of large-scale jets 
in other galaxies, it was not expected that the Milky 
Way's SMBH would produce such a relativistic colli- 
mated structure, given its current quiescence. However, 
the MW must have undergone phases of nuclear activ- 
ity in the past in order for the SMBH to grow, and it 
is plausible that signs of past activity are still visible. 
One might expect relics of past activity in high energy 
cosmic rays (CRs) and hot gas, perhaps far off the disk. 
The sensitivity and angular resolution of the Large Area 
Telescope (LAT) on board the Fermi Gamma-ray Space 
Telescope ( Atwood et al.l l2009) make possible the search 
for inverse Compton (IC) gamma rays from a Galactic 
jet. 

In this work, we use LAT data at 0.3 GeV < E 1 < 100 
GeV to look for unexpected diffuse Galactic gamma-ray 
structure. We will show evidence for a large-scale col- 
limated double-jet structure, which appears symmetric 
with respect to the Galactic center (GC). For the re- 
mainder of this paper, we refer to these collimated jet- 
like structures as the "gamma-ray jets" or simply the 
"jets". We also argue that the southern jet has pro- 
duced a large cocoon feature, visible over a wide range 
of gamma-ray energies. These jets and cocoon features 
might be associate d with the prev iously discovered Fermi 
bubble structures (|Su et al.ll2010[ h and may hold the key 
to understanding their origin. 

In Section [21 we describe the Fermi-L AT data selection 
and our data analysis procedure including map making. 
In Section [3l we describe our model of the Galactic dif- 
fuse gamma-ray emission. In Section [H we show that 
the gamma-ray maps constructed from three-year Fermi- 
LAT observations reveal evidence of large-scale Galactic 
jet features along with a south cocoon structure. We 
characterize the morphology of the jet/cocoon system 
in some detail and employ regression template fitting to 
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determine the jet and cocoon spectra in Section [5l We 
calculate and discuss the radio luminosity of the jet in 
Section [6] and carefully discuss the statistical significance 
of the jet structure in Section [7| Finally, we discuss the 
implications of the presence of the Galactic jet and future 
observations in Section [U 

2. MAP CONSTRUCTION FROM FERMI-LAT DATA 
2.1. Fermi Data Selection 
The LAT ([Gehrels fc Michelsonl Il999l : lAtwood et al 



significantly lower backgro und compared t o the previ- 
ously released v2_3 maps in Su et al.l ( 2010l ). The "three 
year maps" exclude some short time periods, primarily 
while Fermi passes through the South Atlantic Anomaly. 
We construct maps using the newly released photon 
event list corresponding to the P7_V6 Instrument Re- 
sponse Fu nctions (I RFs) 



20Q9h is the principal scientific instrument on board the 
Fermi Gamma-ray Space Telescope. The point-spread 
function (PSF) is about 0.8° for 68% containment at 1 
GeV and decreases with energy as r^s ~ £? -0 ' 8 , asymp- 
toting to ~ 0.2° at high energy. It is designed to survey 
the gamma-ray sky in the energy range from about 20 
MeV to several hundreds of GeV. 

Using standard photon event selection (zenith angle 
< 100°, etc.) we generate full-sky maps of counts and 
exposure using HEALPix, a convenient equal-area iso- 
latitude full-sky pixelization widely used in the cosmic 
microwave background (CMB) community. 4 Spherical 
harmonic smoothing is straightforward in this pixeliza- 
tion, and we smooth each map by the appropriate ker- 
nel to obtain a Gaussian PSF of 2° FWHM. Because 
the PSF of the initial map must be smaller than this, 
at energies from 300 MeV to 1 GeV we use only front- 
converting events, which have a smaller PSF. Above 1 
GeV, we combine front- and back-converting events. 

We make use of the Fermi Pass 7 (P7) data products, 
using the latest publicly available event reconstruction 
algorithms 5 . For the purpose of studying Galactic diffuse 
gamma-ray emission with minimal non-photon contam- 
ination (low background), we select events designated 
ULTRACLEAN 6 class that have the most stringent data se- 
lection criteria, so that any diffuse features appearing 
on the gamma-ray maps are not due to CR contamina- 
tion. Photons coming from the bright limb at Earth's 
horizon, dominantly produced by grazing-incidence CR 
showers coming directly toward the LAT are a strong 
source of contamination. We minimize this background 
by restricting to events with zenith angle less than 100° 
as suggested in the Fermi Cicerone 7 . 

2.2. Map Making 

Our current gamma-ray maps (v3_3) constructed from 
the three year Fermi data have greater signal/noise and 

4 HEALPix software and documentation can be found at 
http://healpix.jpl.nasa.gov, and the IDL routines used in 
this analysis are available as part of the IDLUTILS product at 
http : / / sdss3data . lbl . gov/ software/ idlutils. 

5 Details at http://fermi.gsfc.nasa.gov/ssc/ 

6 Maps produced with CLEAN have about 20% more events. Re- 
peating our analysis on these maps yields similar results, with sim- 
ilar error bars at E < 10 GeV, and slightly larger error bars at 
E > 10 GeV, where background is relatively more important. 

7 http: //fermi .gsfc .nasa.gov/ ssc/dat a/ analysis/document at ion/ 



As 



Su et al. 



( 2010[ ), we construct maps of 
front-converting and back-converting events separately, 
smooth to a common PSF, and then combine them. To 
reveal the diffuse emission, we subtract point sources 
using the Second Fermi-L AT catalog (2FGL), which is 
based on 24 months of LAT observations 8 and the P7_V6 
event selections and IRFs. The PSF and effective area of 
the Fermi-LAT vary with energy, and we subtract each 
point source from the maps in each energy bin, using the 
in-flight version of the PSF contained in the P7_V6 IRFs. 
We produce the exposure maps using the gtexpcube task 
in the Fermi Science Tools. 

In principle the PSF of each event depends on the an- 
gle between its arrival direction and the instrument axis. 
The LAT scan strategy results in a small variation of the 
PSF with position on the sky. We use the in-flight deter- 
mination of the PSF that neglects these details 9 . For the 
400 brightest and 400 most variable sources, the subtrac- 
tion is noticeably imperfect, so we interpolate over the 
core of the PSF after subtracting the best estimate. We 
take care to expand the mask for very bright sources 
(Geminga, 3C 454.3, and LAT PSR J1836+5925) and 
large sources (Orion and the Magellanic Clouds). The 
resulting map is appropriate for the analysis of diffuse 
emission at \b\ > 3°. At |6| < 3° the maps are severely 
compromised by the poor subtraction and interpolation 
over a large number of point source s. Further de t ails o f 



the map proc e ssing may be found in iDobler et al.l (|2010h 



and ISu et al.l ([20101 ). The v3_3 maps used in this work 
are available for download 10 in both FITS and jpeg for- 
mats. 

3. DIFFUSE GALACTIC GAMMA-RAY EMISSION MODELS 

At low (~1 GeV) energies, and near the Galactic plane 
(\b\ < 20°), the gamma rays observed by Fermi-LAT are 
dominated by photons from the decay of 7r° particles, 
produced by the collisions of CR protons with ambient 
gas and dust in the interstellar medium (ISM). Colli- 
sions of CR electrons with the ISM (primarily protons, 
but also heavier nuclei) produce bremsstrahlung radia- 
tion. The CR electrons also IC scatter the interstellar 
radiation field (ISRF) up to gamma-ray energies. In 
order to reveal the gamma-ray jet features, significant 
7T° emission, bremsstrahlung, and IC emission from the 

8 http: //fermi .gsfc .nasa.gov/ssc/data/access/lat/2yr_catalog, 
the file we used is gll_psc_v05 . f it 

9 See http : / /fermi . gsfc . nasa . gov/ ssc/ data/ analysis/ documentation 
/Cicerone/Cicerone_LAT_IRFs/IRF_PSF . html 

10 Available at http : / /fermi . skymaps . info 
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Fig. 1. — A large scale gamma-ray cocoon feature is revealed by template decomposition of the gamma-ray maps in energy range of 1 — 2 
GeV (upper-left), 2 — 5 GeV (upper-right), 5 — 10 GeV (lower-left), and 10 — 50 GeV (lower-right) maps constructed from Fermi-L AT 
three year observations. Bright point sources have been masked and fainter ones subtracted. The maps are smoothed by a Gaussian kernel 
with FWHM of 90 arcmin. We remove the dominant diffuse Galactic signal from 7r° (and bremsstrahlung) gam ma-rays prod uced by CR 
protons (electrons) interacting with the ISM using a dust map template constructed based on far IR data (Schleg el et al.lll998h . A smooth 
disk model is also subtracted to reveal the structure deeper into the plane. This model mostly removes the IC gamma rays produced by 
CR electrons interacting with the interstellar radiation field including CMB, infrared and optical photons. We also subtract a uniform 
template of the Fermi bubbles to remove the corresponding gamma-ray em ission. The cocoon feature is revealed on the east (left) side of 
the previously discovered southern Fermi bubble structure (|Su et al.ll2010h . with relatively sharp edges. 



Galactic disk must be removed. We take two approaches 
to this foreground removal: one is to use the Fermi Dif- 
fuse Galactic Model 11 provided by the Fermi team; the 
second is to employ a linear combination of templates of 
known emission mechanisms, using existing maps from 
multi- wavelength observations and/or constructed geo- 
metric templates. 

3.1. Fermi Diffuse Galactic Model 

The LAT Diffuse Galactic Model is a comprehensive 
model of Galactic gamma-ray emission from the ISM 12 . 
The Fermi diffuse model is primarily designed as a 
background template for point source analysis or in- 
vestigation of small-scale diffuse structures, and comes 
with a number of caveats. However these caveats ap- 
ply mainly near the Galactic plane, and at E > 50 GeV. 

11 http: //fermi .gsfc .nasa.gov/ssc/data/ 

12 Available from fermi.gsfc.nasa.gov/ssc/data/ 



It is nevertheless useful for qualitatively revealing fea- 
tures in the diffuse emission at high latitude. In this 
work, we use the version of Fermi diffuse Galactic model 
gal_2yearp7v6_v0 .fits. 

3.2. Simple Template- based Diffuse Galactic Model 

The 7T° /bremsstrahlung gamma-ray intensity is pro- 
portional to the ISM density x the CR proton/electron 
density, integrated along the line of sight. As long as 
the CR proton/electron spectrum and density are ap- 
proximately spatially uniform, the ISM column density 
is a good tracer of the resulting gamma-ray distribution 
from 7T° /bremsstrahlung emission. Because the emission 
is optically thin, subtracting a simple template of the 
7T° gammas helps to reveal the gamma-ray jets/cocoon 
along with the Fermi bubbles, especially towards the 
Galactic plane. The ISM column density is expected 
to be strongly correlated with other tracers of the ISM. 
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Our previous work ( Su et al. 2010h used the SFD map 



of G alactic dust, based on far IR data ([Schlegel et al 
1998h . The dust map has some advantages over gas 



Fig. 2. — We show gamma-ray maps in Galactic coordinate (—50° < b < 50°, —25° < t < 25°) from Fermi-LAT for photon energies 
0.8 — 3.2 GeV. The Galactic center is marked with a cross sign in the center of the maps. Two large-scale collimated jet-like features (the 
gamma-ray "jets") are revealed. The right panel shows the same image as the left panel, but with a dashed line representin g the direction 
of the suspected jet. Point sources have been subtracted based on the Second Fermi-LAT catalog (2FGL) (|The Fermi-LAT Collaboration! 
l201lh . and large sources, including the inner disk (—5° < b < 7°), have been masked. The maps are smoothed by a Gaussian kernel with 
FWHM of 2°. The Fermi diffuse Galactic model (Pass7_V6) has been subtracted to remove known large scale diffuse gamma-ray emission. 

We use the SFD dust map as a template of the 
7T° gamma foreground. The correlation between Fermi 
gamma-ray maps and the SFD dust map is striking, and 
most of the ISM emission is removed by this subtraction. 
To reveal the structure de eper into the p lane, a simple 
disk model is subtracted (|Su et al.l l2010l ) 13 . The disk 
model mostly removes the IC gamma rays produced by 
cosmic ray electrons interacting with the ISRF including 
CMB, infrared, and optical photons; such electrons are 
believed to be mostly injected in the Galactic disk by 
supernova shock acceleration before diffusing outward. 
The two jets are symmetric with respect to the Galactic 
plane and well centered on the GC, thus they are unlikely 
to be local structures. 



maps. One shortcoming of using Hi and CO maps is 
the existence of "dark gas" (|Grenier et al.ll2005f ). clouds 
with gamma-ray emission that do not appear in the Hi 
and CO surveys. Th ese features are seen in dust maps 
( Schlegel et al.l 1 19981 ) and may simply be molecular H 
clouds under- abundant in CO. For the dust map, there 
are no problems with self absorption, thus no concerns 
about "dark gas," and the SFD dust map has sufficient 
spatial resolution (SFD has a FWHM of 6', and LAB 
Hi is 36') . On the other hand, SFD contains no velocity 
information, so it is impossible to break the map into 
Galactocentric rings as is done with Hi and CO maps. 
Nevertheless, it is instructive to employ the SFD map 
to build a very simple Galactic diffuse gamma-ray tem- 
plate. The goal is to reveal the gamma-ray jet/cocoon 
structures by removing the expected diffuse emission in 
a fashion that makes as few physical assumptions as pos- 
sible. 



4. GAMMA-RAY COCOON AND EVIDENCE OF GALACTIC 

JET 

The gamma-ray cocoon in the south is visible in the 3- 
year LAT maps at energies from 1 to 50 GeV (Figure [Tj). 



13 The functional form of this disk template is (esc |6|) 
latitude and a Gaussian {pn = 30°) in longitude 
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Fig. 3. — Left panel: spectral energy distribution of the components in our six-template fit. The dust-correlated spectrum (red 
short-dashed line), traces 7r° emission. The disk-correlated emission (green dashed), approximately traces the soft inverse Compton and 
bremsstrahlung components. The spectrum of the uniform emission (dotted brown line) includes the isotropic part of the extragalactic 
background and cosmic-ray contamination. The spectrum of emission correlated with Loop I (solid orange) has a spectrum similar to the 
disk-correlated emission. In contrast to these soft-spectrum components, the Fermi bubble template (blue dot-dashed) and the gamma-ray 
cocoon (black solid) have notably harder (consistent with flat) SEDs. Vertical bars show the marginalized 68% confidence range derived 
from the parameter covariance matrix for the template coefficients in each energy bin. Right panel: same, but with one additional template 
representing the gamma-ray jet. For clarity, we only show the spectrum of the gamma-ray jet, the cocoon and the Fermi bubbles, and 
compare their spectra with the softer 7r° emission (red dashed). 

By inspection, the major symmetric axis of the cocoon 
lines up with the Galactic center to within ~ 5°, and the 
major to minor axis ratio is about a factor of three. The 
gamma-ray luminosity of the south cocoon with Galactic 
latitude \b\ > 20° and in the energy range 1 — 100 GeV is 
(5.5±0.45) x 10 35 erg/s. The morphology of the gamma- 
ray cocoon resembles the observed ra dio cocoon struc- 
tures of Fanar off- Riley (FR) type II ([Fanaroif fc Rilev 
1974h active radio galaxies (e.g. Cygnus A), which have 
been found surrounding collimated large-scale jets and 
may be formed by backflow of magnetized jet plasma. 
How the gamma-ray cocoon system formed and retains 
its tight columnar shape despite traveling from the GC 
for ~10 kpc is intriguing. The presence of this large scale 
cocoon suggests collimated injection of high energy par- 
ticles from the inner Galaxy. 

In Figure [21 we show the 0.8 — 3.2 GeV gamma-ray map 
of the inner Galaxy. A pair of collimated linear features is 
revealed in Figure [21 with similar morphology in each en- 
ergy bin. There are no other apparent large scale features 
in the residual maps. The gamma-ray jets do not appear 
to be associated w ith the well-known Loop I structure 
()Ber khui j sen! 1 1 9 73[ ) . The north jet extends from the GC 
to Galactic coordinates (£,b) = (—11,40), and the south 
jet extends from the GC down to (£,b) = (11, -44). Al- 
though the jets are faint, three lines of evidence suggest 
that they are real: (1) the jets both emanate from the 
GC, in nearly opposite directions; (2) the jets extend 



away from the GC to about the ed ge of the prev iously 
discovered Fermi bubble structure flSu et al.ll2010h : and 
(3) the south jet aligns with the symmetry axis of the 
cocoon structure. The morphology of the jet and cocoon 
gamma-ray feature and the possible association with the 
Fermi bubbles strongly suggests the Galactic center ori- 
gin of these structures and implies recent activity toward 
the inner Galaxy. 

On a much smaller scale, Chandra X-ray observations 
have revealed a faint jet -like feature at su bparsec scale, 
pointing toward Sgr A* (|Muno et al.l 1200 8). It is tempt- 
ing to interpret this feature as a jet of synchrotron- 
emitting particles ejected from the SMBH. However, due 
to it's much smaller size and ~15° misalignment with 
the gamma-ray jet, there is no clear association. Lim- 
its on the gamma-ray jet from other wavelengths do not 
strongly constrain its nature, underscoring the need for 
further mult i- wavelength studies. 

The projected direction of the gamma-ray jets is about 
15° from the north-south axis of the Galaxy. If the 
gamma-ray jets constitute the projection of double large- 
scale jets symmetric to the Galactic plane and the GC, 
taking the distance from the solar system to the GC Rq 
= 8.5 kpc, each projected jet is ~10 kpc. The small angle 
of the projected jet relative to the north-south Galactic 
axis suggests that the apparent length of ~10 kpc is a 
good approximation of the spatial scale of the Galactic 
jet. The width of the jets appears to be < 5°. The to- 
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Fig. 4. — This figure shows the best fit linear combination maps corresponding to the spectra in Figure[T] To increase signal/noise, larger 
energy bins are used. The left column shows the linear combination of the disk, Loop I, uniform, bubble template, the south cocoon, and the 
south jet template, that provide the best fit to the Fermi maps (middle column) after subtracting the best fit SFD dust template. Because 
the 7T° emission traced by SFD is so bright, it is subtracted from both the models and data shown in this figure. The difference maps (data 
minus template model) are in the right column. The template fitting is done for the region with \b\ > 20° to avoid contamination from 
the Galactic disk. The subtraction of the model largely removes the features seen in the Fermi maps with \b\ > 20°. We have also masked 
the inner Galactic plane region (|6| > 4° and |Z| > 60°), which is significantly contaminated by point sources. We use the same gray scale 
for all panels. We find that both the disk IC template and Loop I features fade away with increasing energy, but the jet template does 
not. The oversubtraction in the residual maps, especially in the lower energy bins, is due to the simple disk IC model, which is not a good 
template across the entire disk. However, in the fit region (|6| > 20°), the residual maps are consistent with Poisson noise without obvious 
large scale features. 
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Fig. 5. — This figure is the same as Figure [3] except with the 
north jet template added to the regression. The north and south jet 
spectra are noisy, but consistent with each other, and consistently 
hard. Note that the upper limit is shown for 3cr. The spectra 
plotted here are available in Tables 1 and 2. 

tal luminosity of the north and south jet-like features is 
(1.8 ± 0.35) x 10 35 erg/s at 1 - 100 GeV (see Table 2). 
The jets do not align with any plausible artifacts relat- 
ing to the Fermi orbit or scan direction, or any other 
known systematics. The north and south jets have sim- 
ilar integrated gamma-ray flux, so there is no evidence 
that the jets are close to the line of sight; otherwise, 
the approaching jet would appear substantially Doppler 
brightened relative to the receding one. 
This gamma-ray jet/cocoon structure is not visible in 



the R O SAT All-Sky Survey of soft X-ray (jSnowden et al 



19971 ). Wilkinson Microwave Anisotropy P robe (WMAP) 



microw ave maps, radio maps at 408 MHz (Ha slam et al 



1982ah . or any other available radio maps. We estimate 
expected radio spectra in Figure [9] corresponding to elec- 
tron spectral indices 1.5 < 7 < 2.5, and find that over 
this range, current full-sky radio surveys are not sensi- 
tive enough to detect the jet. However, combining the 
gamma-ray energy spectrum along with the radio limits, 
we find that the CR electrons producing the jet have a 
hard spectrum with spectral index > — 2. 

We now turn to the relation between the jet/cocoon 
and the Fermi bubbles. We recently discovered two 
giant gamma-ray bubbles with a total luminosity from 
1 to 100 GeV of 2.0 x 10 37 erg/s 14 , extending - 50° 
above and below the GC, with a width of ~ 40° in lon- 
gitude, and found them to be spatially correlated with a 
hard-s pectrum microwav e excess (known as the WMAP 
haze) (jFinkbeinerl l2004f ) and large-scale X-ray features 



14 We note that in lSu et all (|2010h , the total gamma-ray lumi- 
nosity of the Fermi bubbles was mis-quoted as 4.0 x 10 37 erg/s. 



( Su et alllioioh . Galactic shock waves produced by en- 
ergetic explosions at the Galactic nucleus or by a high 
rate of supernova explosions in the nuclear disk will be 
channeled by the decreasing ambient gas density in di- 
rections perpendicular to the Galactic plane. 

The gamma-ray emission associated with these bub- 
bles has a significantly harder spectrum (dN/dE ~ E~ 2 ) 
than the IC emission from electrons in the Galactic disk, 
or the gamma-rays produced by decay of 7r° from proton- 
ISM collisions. The Fermi bubbles are likely formed dur- 
ing an active phase in the Galactic center ~ 10 6 — 10 7 
years ago with jets ejected from the central supermas- 
sive black hole. The bubble region might consist of de- 
celerated jet material, radiating isotropically. Observa- 
tional data and numerical simulations indicate that the 
energy required to form the bubbles is on the order of 
10 5 



55 1Q 58 erg |j Su et al J I2010I : iGuo fc Mathewsl 12011 



Guo et aTll2011 ). Even though the current jet luminosity 
is relatively faint, the discovery of the jet/cocoon system 
generally supports the AGN hypothesis for the origin of 
the bubbles. 



5. ENERGY SPECTRUM OF THE GAMMA-RAY COCOON 
AND JETS 

We have shown in Figure [T] that the gamma-ray jets, 
the cocoon, and the Fermi bubbles all have a harder 
spectrum than other large-scale diffuse ISM emission. In 
order to measure the hardness of the spectrum and to ex- 
plore possible mechanisms responsible for the observed 
gamma-ray emission, we maximize the likelihood that 
the maps are described by a linear combination of spa- 
tial templates, one for each gamma-ray emission mech- 
anism. We fit a coefficient to each emission component 
using a multi-linear regression of simple templates, one 
energy bin at a time. By combining results from 12 
logarithmically-spaced energy bands from 300 MeV to 
300 GeV, we determine a spectral energy distribution for 
each component. In each fit, we maximize the Poisson 
likelihood of a simple diffuse emission model involving 
the seven following templates : the SFD map of Galac- 
tic dust ([Schlegel et al.l fl99 8). the simple disk model, 
the Fermi bub bles, a template of the Loop I structure 
()Su et al.ll2010f h the gamma-ray cocoon, the gamma-ray 
south jet, and a uniform background as template to ac- 
count for backgrou nd gamma-ray emission and cosmic 
ray contamination ( Su et al. 2010l ). All maps and tem- 
plates have been smoothed by a Gaussian kernel with 
FWHM of 2° for the regression analysis. Systematic un- 
certainties are dominated by the imperfect representation 
of the diffuse emission by these simple templates. We feel 
that this analysis is sufficient for a rough characterization 
of the cocoon and jet spectra, and provides motivation 
for more thorough analysis using a more physical model 
in the future. 

For each set of model parameters, we compute the Pois- 
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00.5-01.7GeV 




Fig. 6. — Same as Figure [4] except we include the north jet template. 



son log likelihood, 

In C = kj In ji j- m - \n(ki\), (1) 

i 

where fii is the synthetic counts map (i.e., linear combi- 
nation of templates times exposure map) at pixel z, and 
k is the map of observed counts. The last term is a func- 
tion only of the observed maps. We compute errors in 
the Gaussian approximation by inverting the matrix of 



second partial derivatives of — In C to obtain the covari- 
ance matrix, and taking the square root of the diagonals. 
The la Gaussian error corresponds to A In £ = 1/2. We 
refer to Appendix B of Dobler et al. ( 2Q10h for more de- 
tails of the likelihood analysis. Maps of the models are 
constructed from linear combinations of these templates, 
and the (data minus model) residual maps in various en- 
ergy bins are shown in Figured! In this fit, we mask out 
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Fig. 7. — Same as Figure [T] except we add the north jet template 
and split the Fermi bubble into north and south bubbles separately. 
The south cocoon template is included in the fit. The north and 
south jet spectra are consistently hard. The fitting results are not 
significantly affected by the bubble template splitting. 

all pixels with Galactic latitude \b\ < 20°. 

So far we have described the procedure for obtaining 
the correlation coefficient for each template at each en- 
ergy. However, some templates have units (e.g., the SFD 
dust map is in magnitudes of Eb-v reddening) so the 
correlation coefficient has unusual units (e.g. gamma- 
ray emission per magnitude). In such a case we multi- 
ply the correlation spectrum by the average SFD value 
in the (\b\ > 20°) bubble region, to yield the average 
spectrum of dust-correlated gamma-ray emission in this 
region. For the uniform, Loop /, cocoon, jet, and bubble 
templates, no renormalization is done. These templates 
are simply ones and zeros (smoothed to the appropriate 
PSF), so the south jet spectrum is simply the spectrum 
of the south jet template shown in e.g. left panels of 
Figure [4j not the mean of this template over the whole 
bubble region. 

The energy spectra for 7r° emission, bremsstrahlung 
and IC scattering can be calculated using a sample 
GALPR0P model 15 (tuned to match locally measured pro- 
tons and anti-protons as well as locally measured elec- 
trons at ~ 20 — 30 GeV), as an indication o f the ex- 
pecte d spectral shapes. We have shown in (ISu et al 



2010l ) that the energy spectra for the SFD and the sim- 
ple disk template reasonably match the model expecta- 
tions. The dust map mostly traces the 7r° emission, and 
the simple disk model resembles a combination of IC 

15 GALPR0P is a cosmic ray propagation code. It calculates the 
steady state solution to the diffusion-energy-loss equation, given 
the 3D gas distribution, interstellar radiation field, B-field model, 
CR diffusion assumptions, and many other input parameters 
([Strong fe Moskalenko 1999). See http: //galprop. Stanford. edu 



and bremsstrahlung emission. The spectrum for emis- 
sion correlated with the gamma-ray jet and cocoon is 
clearly significantly harder than either of these compo- 
nents, consistent with aflat spectrum in E 2 dN/dE. This 
fact coupled with the distinct spatial morphology of the 
jet/cocoon system indicates that if these gamma-ray fea- 
tures are generated by IC scattering of interstellar ra- 
diation fields, then a separate electron population must 
exist in the jet/cocoon. We also note that the spectrum 
of the jet/cocoon template does not fall off significantly 
at energy <1 GeV as the bubble spectrum does. The 
fitting coefficients and corresponding errors of each tem- 
plate are listed in Table 1 and 2. 

The hardness of this spectrum may be the key to 
deci phering the o rigin of the gamma-ray jet/cocoon. 
In (|Su et alJlioioh . we show that at lower energy (E < 
1 GeV) the bubble spectrum falls sharply with decreas- 
ing energy (becomes dramatically harder than —2). With 
reduced statistical error from three year LAT data and 
the new Pass 7 ULTRACLEAN event selection, we confirm 
significant falling of the bubble spectrum at lower en- 
ergy but not the jet and cocoon components. This is 
true whether we perform the entire analysis with front- 
converting events only, or use front-converting at E < 1 
GeV and both front- and back-converting at E > 1 GeV 
as usual. 

The null hypothesis of zero intensity of the north 
(south) jet is ruled out by 3.1a (4.1a), respectively, and 
5.2a jointly for the whole jet structure. The same fit 
simultaneously finds the cocoon with 12a significance. 
Since the region south of the Galactic center has less 
foreground emission than the north (where p Oph is), we 
focus on the south jet/cocoon for our analysis. Figure [3] 
(right panel) shows the resulting energy spectrum of the 
gamma-rays associated with the cocoon and the south 
jet. Like the Fermi bubbles, the energy spectrum of the 
jet /cocoon is harder than other diffuse gamma-ray com- 
ponents, although the cocoon has a spectrum at < 1 GeV 
different from the Fermi bubbles, which suggests its ori- 
gin from a distinct population of electron CRs. The co- 
coon spectrum is consistent with E 2 dN/dE ~ constant, 
and the north and south jet structures have an energy 
spectrum of E 2 dN/dE ~ E°- 2±0 - 2 at latitudes \b\ > 20°. 
The correlation coefficients for the SFD map and simple 
disk model are multiplied by the average value of these 
maps in the bubble region with a \b\ > 20° cut to obtain 
the associated gamma-ray emission. Given how hard the 
cocoon/jet spectra are up to ~ 100 GeV, bremsstrahlung 
can be ruled o ut as the emission mechanism based on the 
arguments in ( Kino et al. 2QQ9h . 

The gamma-ray jets appear to possess north-south 
symmetry in morphology without noticeable difference 
in intensity. To investigate whether there is any spectral 
difference between the north and south jets, we augment 
the 7-template fit (dust, simple disk, Fermi bubble, co- 
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Fig. 8. — Same as Figure [T] except we add the north jet template and split the bubble template to north and south separately. Note 
that the south cocoon template is included in the fitting. The north and south jet spectrum are consistently hard. 



coon, south jet, uniform, Loop I), with a north jet tem- 
plate. We show the resulting spectrum in Figure [5] and 
the maps in Figure [6l We then repeat the previous fit- 
ting procedure involving the simple disk IC template, but 
splitting the bubble template into north and south bub- 
ble templates, and allow an independent fit of the two 
jet templates along with the two bubbles. The goal is 
to identify variations in the intensity and spectral index 



between the northern and southern jets and the depen- 
dence of the jet spectrum on the combination of different 
emission templates. Even given this freedom, no signif- 
icant spectral differences are found between the north 
and south jets or the north and south bubbles. There 
was no significant improvement of the likelihood for such 
splitting. The spectrum and the maps are shown in Fig- 
ure [3 and Figure [8j respectively. Our conclusion is that 
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the gamma-ray jets appear to be north-south symmetric, 
both with a hard spectrum. This statement is largely 
independent of our choice of template for the disk IC 
emission. In summary, we found no apparent difference 
between the north and the south jets both in morphol- 
ogy and in energy spectrum, which indicates they might 
share the same origin. 
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Fig. 9. — The estimated spectrum of inverse Compton gamma- 
rays (upper panel) and corresponding synchrotron radiation (lower 
panel) originating from a hard electron spectrum along a line of 
sight 4 kpc below the Galactic plane (i.e. b w —25°). The steady- 
state electron spectrum is taken to be a power law, dN/dE oc , 
with index 7 = 1.5 (solid black), 1.7 (blue dashed), 2.0 (green dot- 
ted), and 2.5 (red dash-dotted) in both the upper and lower pan- 
els. In all cases the cosmic ray electron has a range of [0.1, 1000] 
GeV. The interstellar radiation field model is taken from GALPR0P 
version 50p, and the magnetic field is assumed to be 5 fiG for syn- 
chrotron calculation. The data points in the upper panels show 
the south cocoon emission the same as in Figure [3] The arrows 
shows 3cr upper limits rather than data points with la error bars, 
due to the large uncertainties in those energies. The data point 
in the lower panel shows the magnitude of the WMAP haze aver- 
aged over b = —20° to —30° in the 23 GHz K-band, and the gray 
area indicates the ra nge of synchrotron spectra l indices allowed 
for the WMAP haze (|Dobler fe Finkbeinerll2008l ). Two arrows at 
lower frequency shows the 3cr upper limits for the cocoon radio 
emission from Rhodes/HartRAO 2.326 GHz radio c ontinuum sur- 
vey (I Jonas et al .11 19981 ) and Haslam 0.408 GHz map (|Haslam et al.l 
Il982bh . The same population of cosmic ray electrons can consis- 
tently both the radio/microwave observations and the gamma-ray 
observations. 



6. EXPECTED RADIO LUMINOSITY OF THE JET 



In order to estimate the expected radio luminosity of 
the jet, we must either have a physical model of the ra- 
dio and gamma-ray emission mechanisms, or a proto- 
type object for comparison. Very Long Baseline Array 5 
GHz observ ations of blazars dete cted by LAT provide one 
benchmark (|Linford et al.l 120 111 ). This sample of blazars 
shows a median 5 GHz flux of ~ 300 mJy and a median 
gamma-ray flux (0.1 — 100 GeV) of 3 x 10 -8 photons 
cm _2 s _1 , although individual objects can deviate from 
this ratio by an order of magnitude or more. BL Lac ob- 
jects have a spectrum of dN/dE ~ E~ 1,8 in this energy 
range, and therefore have 1/15 of their 0.1 < E < 100 
GeV photons in the range 1 — 2 GeV. Assuming such 
a spectrum, a 5 GHz flux of 300 mJy corresponds to 
4 x 10- 9 GeV cm^s- 1 (1-2 GeV). Therefore, the 5 
GHz radio signal corresponding to our observed 4 x 10 _T 
GeV cm -2 s -1 sr -1 would be 30 Jy/sr, much smaller than 
our radio limits in Figure However, because BL Lac 
objects are thought to contain a relativistic jet pointed 
along the line of sight, they may be a poor analog for the 
Galactic jet. 

An alternative is to take a magnetic field constrained 
by radio and F araday rotation measure observations 
( Sun et aT]l2008 ) and an interstellar radiation field (taken 
from GALPR0P, 4 kpc below the Galactic center) to com- 
pute the IC gammas and synchrotron given a power-law 
CR electron spectrum. The high-energy cutoff of the 
electron spectrum is important for the IC gammas, es- 
pecially for spectra harder than dN/dE ~ E~ 2 . As an 
example, we have computed the IC gamma-ray and syn- 
chrotron emission for 4 power laws, with spectral index in 
the range 7 = 1.5 — 2.5 with a cutoff at ^electron = 1 TeV 
(Figure [9]). The electron spectra are normalized such 
that the IC gammas go through the data points. The 
synchrotron is then compute d assuming a 5 11G magnetic 
field and the GALPR0P ISRF ([Porter et al.ll2008l ) at 4 kpc 
below the Galactic p l ane. U pper limits from maps at 408 



MHz (Haslam et al 



1982bh 



and 2.3 GHz (Jonas et al 



19981 ). and a data point from WMAP 23 GHz 4 kpc off 



the disk (whic h together serve as an upp er limit for the 
jet brightness (|Dobler fc Finkbeinerlboolj )) are included. 
Most of the lines pass below the limits, making it plausi- 
ble that the radio signal corresponding to a gamma-ray 
jet would not yet have been observed. Additional uncer- 
tainty arises from the fact that the 1 GHz radio signal is 
produced by order 1 GeV electrons, whereas the gamma- 
ray signal arises from order 100 — 1000 GeV electrons. If 
the electron spectrum has any downward curvature over 
this range in energy, the radio /microwave limits become 
even less constraining. Another uncertainty comes from 
the choice of energy cutoff for the power law. We con- 
sidered spectra with cutoffs at 3 and 5 TeV with index 
7 ~ 3 to fit the gamma-ray data. Such spectra predict 
a synchrotron signal bright enough to be ruled out by 
the radio data. This implies that the required CR elec- 
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Fig. 10. — Full sky residual maps after subtracting the Fermi diffuse emission model from the Fermi-L AT three year gamma-ray maps. 
Upper panel: point sources have been subtracted, and large sources, including the inner Galactic disk, have been masked. The Fermi 
bubble structure has been included in the diffuse emission model. We find no large-scale jet features other than the central Galactic jet 
toward the inner Galaxy. The Galactic jet feature is crossing the Galactic center, and aligns with the long axis of the cocoon structure. 
Lower panel: the Fermi bubble edge has been marked in green dashed circles above and below the Galactic center, overplotted on the 
gamma-ray map the same as the upper panel. Both the north and the south jets approximately end on the edge of the Fermi bubble 
structure. 



tron spectrum is much harder than the locally measured 
index of 7 ~ 3 — 3.25. 

7. STATISTICAL SIGNIFICANCE OF THE JET 

In Figure \W\ we show the full sky residual maps at 
0.8—3.2 GeV after subtracting the Fermi diffuse emission 
model. This figure is the full sky version of Figure [H We 
identify a pair of Galactic jet features toward the inner 
Galaxy, present both above and below the Galactic plane 
and aligned with the Galactic center. There are no other 
significant large-scale jet structures appearing in the full 
sky gamma-ray map. 

Our best-fit values are derived from a maximum like- 
lihood analysis, in which we maximize the Poisson prob- 



ability of observing the observed counts, given a set 
of model parameters. The error bars are derived from 
the parameter covariance matrix. In our previous work 
( Su et all l2oTol ). we analyzed numerous mock maps to 
verify that our code gives unbiased results with the cor- 
rect uncertainties, at least in the case where the model 
(plus Poisson noise) is a good description of the data. 
However, our fits contain systematic residuals, and it is 
necessary to investigate to what extent these might be 
able to mimic a signal. To this end, we have repeated 
the 7-component fit with the cocoon and jet structures 
shifted around the sky at 5 deg intervals in longitude 
and determined the best fit coefficient for each case by 
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Fig. 11. — As a verification of our uncertainty estimate for the 
Galactic jet structure, we have rotated the south jet and cocoon 
templates in longitude around the sky, and determined the best 
fit coefficient for each case. This histogram shows the distribution 
of the coefficients overplotted with the fitted Gaussian distribu- 
tion. The coefficient distribution is well centered around zero and 
well fitted with a Gaussian noise distribution. The south jet tem- 
plate with a coefficient 0.182 is about Act significance from the 
background noise estimation. The estimated uncertainty of the jet 
template agrees well with our estimation of the significance of the 
south jet based on the Poisson likelihood fitting. 

fitting the same templates as in Figure El By placing the 
templates at different position along the longitudes where 
we expect no significant large scale diffuse structures, we 
can measure the rms of the jet and cocoon coefficients, 
and this rms may exceed that due to the Poisson noise. 
In some parts of the sky, (e.g. £ ~ 290 — 320°), the 
7T° emission is somewhat oversubtracted, which pushes 
the cocoon coefficient negative, effectively renormalizing 
the sky background for the jet. As expected, the distri- 
bution of jet template coefficients (shown in Figure ITTj) 
is centered around zero, consistent with a Gaussian of 
mean zero and a = 0.044. The Galactic jet coefficient is 
more than 4a out of the background noise distribution, 
which is consistent with the estimated significance of the 
south jet derived from the covariance matrix of our Pois- 
son likelihood fit. We interpret this agreement to mean 
that systematic errors in foreground modeling have not 
substantially distorted the meaning of the south jet's 4a 
formal significance. 

8. DISCUSSION AND CONCLUSION 

Given the alignment of the cocoon with the southern 
bubble edge, the cocoon/jet is probably associated with 
the bubble, and their origins may be intertwined. A pos- 
sible scenario in which a jet creates both the cocoon and 
the Fermi bubble is as follows. The propagation of a 
large-scale relativistic jet may generate a double bow- 



shock structure at t he head of the jet ( Blandford fc Reesl 
19741 : IScheuer1ll974h . Entrained energy and matter are 



pushed aside due to a high pressure gradient and cre- 
ate a hot cocoon around the jet. The cocoon applies 
sufficiently high pressure to collimate the jet and sub- 
stantially reduce its opening angle. Continuous injection 
of relativistic electrons can be produced in the forward, 
reverse, and re-collimation shocks emitting synchrotron 
radiation and gamma-ray IC emission. The jets deliver 
kinetic energy to the surrounding interstellar gas and the 
jet/medium interactions could be strong enough to ac- 
celerate particles and produce non-thermal radiation. In 
this scenario, as the jets emerge from the GC, they ex- 
pand freely until they are recollimated when their ram 
pressure falls to the thermal pressure of the surrounding 
cocoon. The jet decelerates once the accumulated mass 
of the swept up ISM gas becomes similar to that carried 
by the jet. The cocoon and Fermi bubble could be ex- 
plained by jet shocked material (re-confinement region 
and cocoon) and the ambient material could be shocked 
by the bow shock. However, the 15° misalignment be- 
tween the jet/cocoon and the bubbles, and the absence 
of an obvious northern cocoon, suggests that this is not 
the whole story. We cannot rule out the possibility that 
the Fermi bubbles were formed by an earlier energetic 
event, while the jet might be produced more recently 
and be penetrating through the older structure. 

The gamma-ray jet is the first collimated jet struc- 
ture found in gamma-rays and the only collimated jet 
close enough to resolve with Fermi-L AT. The presence 
of such a large-scale Galactic jet in our Milky Way pro- 
vides an ideal nearby laboratory for studying basic ques- 
tions about jet formation, acceleration and collimation. 
It is likely that similar systems are not uncommon. For 
example, the recently observed unusual transient source 
(Swift J164449.3 + 573451) has been understood as a 
newly formed relativistic outflow launched by transient 
accretion onto a SMBH with mass similar to the Galac- 



tic center SMBH ( Zauderer et al. 2011: Burrows et al 



20111 ). It provides the evidence that a normal galaxy 



can transform to an AGN-like phase and produce a rela- 
tivistic jet. Similar structures on much larger scales have 
been found in galaxy clusters with radio lobes inflated by 
the jets of SMBH at the center of powerful radio galaxies 
( McNamara et alJliooih . 

Follow-up observations at other wavelengths are re- 
quired to advance our understanding of the jet and co- 
coon structures, and their relation to the Fermi bubbles. 
In coming years, eRosita will survey the whole sky in 
X-rays with approximately 30 times ROSAT sensitiv- 
ity, in the medium energy X-ray range up to 10 keV 
with an unprecedented spectral and angular resolution 
( Cappelluti et al.l 120111 ). It has the ability to reveal im- 
portant information about the state of the material in 
the jet, cocoon, and bubble structures. Such multi- 
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E Range (GeV) Energy 



Uniform 



SFD Dust Simple Disk 



Bubble Bubble (1.6 yr) South Cocoon 



0.3-0.5 
0.5-0.9 
0.9 - 1.7 
1.7-3.0 
3.0 - 5.3 
5.3-9.5 
9.5 - 16.9 
16.9 - 30.0 
30.0 - 53.3 
53.3 - 94.9 
94.9 - 168.7 
168.7- 300.0 



0.4 
0.7 
1.3 
2.2 
4.0 
7.1 
12.7 
22.5 
40.0 
71.1 
126.5 
225.0 



1.376 ± 0.007 
1.175 ± 0.007 
0.897 ± 0.007 
0.734 ± 0.006 
0.562 ± 0.007 
0.500 ± 0.008 
0.392 ± 0.009 
0.287 ± 0.011 
0.244 ± 0.013 
0.169 ± 0.014 
0.130 ± 0.016 
0.086 ± 0.017 



1.602 
1.696 
1.489 
1.104 
0.778 
0.475 
0.305 



0.019 
0.019 
0.019 
0.016 
0.017 
0.020 
0.023 



0.236 ± 0.026 
0.156 ± 0.030 
0.170 ± 0.036 
0.060 db 0.039 
0.045 ± 0.040 



0.451 
0.458 
0.383 
0.324 
0.249 
0.170 
0.117 
0.111 



0.011 
0.012 
0.012 
0.010 
0.011 
0.013 
0.015 
0.018 



0.028 ± 0.020 
0.014 ± 0.022 
0.014 db 0.025 
-0.010 db 0.027 



0.024 db 0.020 
0.128 ± 0.021 
0.167 db 0.022 
0.290 ± 0.018 
0.295 ± 0.021 
0.363 ± 0.025 
0.365 ± 0.030 
0.307 ± 0.035 
0.279 db 0.041 
0.309 db 0.050 
0.241 ± 0.057 
0.040 ± 0.052 



0.035 
0.211 
0.321 
0.436 
0.353 
0.343 
0.205 
0.263 
0.217 
0.251 
0.319 
-0.015 



0.033 
0.037 
0.044 
0.036 
0.043 
0.049 
0.055 
0.068 
0.083 
0.120 
0.162 
0.194 



0.373 db 0.055 
0.315 ± 0.060 
0.415 db 0.066 
0.264 ± 0.054 
0.321 db 0.064 
0.306 db 0.078 
0.298 db 0.094 
0.473 db 0.119 
0.051 ± 0.125 
0.132 db 0.152 
0.079 db 0.180 
0.071 db 0.173 



TABLE 1 

The template fitting coefficients and errors corresponds to Figure [5] and Figure [6] The gamma-ray luminosity in each 
energy range is shown in the unit of kev cm~ 2 s~ ^r" 1 . for compari son, we also listed the fermi bubble luminosity from 

(ISU ET AL.ll2010h . 



E Range (GeV) 


Energy 


North Jet 


South Jet 


0.3-0.9 


0.6 


0.235 ± 0.139 


0.107 ± 0.112 


0.9-3.0 


1.8 


0.248 ± 0.163 


0.208 ± 0.115 


3.0-9.5 


5.6 


0.248 db 0.170 


0.002 db 0.120 


9.5 - 30.0 


17.6 


-0.036 ± 0.223 


0.318 ± 0.182 


30.0 - 94.9 


55.6 


0.526 db 0.339 


0.656 db 0.266 


94.9 - 300.0 


175.7 


0.151 db 0.446 


0.770 ± 0.395 



TABLE 2 

The template-fitting Coefficients and Errors of the 
North and South jets correspond to Figure [5] and Figure 
\6\ The gamma-ray luminosity in each energy range is 

SHOWN IN THE UNIT OF keV cm~ 2 s~ 1 sr~ 1 . 
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wavelength exploration will advance our understanding 
not only of our own Galaxy, but of black hole accretion 
and jet formation in general. 
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